
*Calculate the own and cross-price elasticity matrix for each option and also each income quartile seperately.

clear all

*read in the data: 
use "$root/Data/Produced/EV_welfare04.dta"


order price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4, after(choice)


sort id car_option

*Set up the matrices and vectors needed for the computation:
mata

ndraws=150

b_f= (0.021183645,-2.593121925,0.154611204,-0.321820168,0.037136405,0.167413791,0.378321531,0.614503527,0.858999103,-0.005445116,-0.010654642,0.136554812,-0.083769393,-0.022510516,0.007285076,0.06246904,1.476449254,2.260850432,0.424973893,0.853075696,0.6579357,1.4180835,0.027265514,0.980066634,0.090419384,0.171157248,-0.142510175,-0.661616614,-0.177476937,0.595728981,0.214934965,0.378106038,-0.080343758,0.394597341,0.734410779,0.078025015,1.019999496,0.291640034,-0.620171333,-0.126470184,-0.114005609,-0.169462474,0.335901857,-0.232550483,0.070071269,0.05306741,0.003621919,0.006019563,0.024778481)'
/*
** Define the estimated vectors of random coefficients:
*/
beta_f=J(1,ndraws,b_f)
rseed(1234)
beta_p = rnormal(1,ndraws,-0.080965469,0.013364931)
beta_he = rnormal(1,ndraws,0.263640717,0.000661958)
beta_we = rnormal(1,ndraws,0.000192522,0.000045627)
beta_var = rnormal(1,ndraws,-0.040212482,0.001878464)

beta= (beta_p \ beta_he \ beta_we \ beta_var \ beta_f)

mata drop beta_he beta_we beta_var b_f


end

*cd "$root\paper\outputs_stata\welfare_calculations04"

** Calculate the initial probability to be able to cross-check the results: 


mata
X=.
st_view(X,.,"price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4")

l_fit=exp(X*beta)



/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id")
V=panelsetup(id, 1)

    ncol   = ndraws
    bigtot = J(rows(id), ncol, .)
    for (i = 1; i <= rows(V); i++) {
        X1 = panelsubmatrix(l_fit, i, V)
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1))
}

mata drop X1 ncol V

_editmissing(l_fit,0)

prob_all=l_fit:/bigtot
mata drop l_fit X 

prob=mean(prob_all')

mata drop prob_all

st_addvar("double", "prob")
st_store(.,"prob",prob')

mata drop prob
end

frame create results_elast
frame create results_elast_w


*Calculate the elasticitities by simulating a price increase for each option by 1%: 
forvalues x=1/489 {
	
	
timer on 1
/* Simulate a 1% price increase:
		gen price_alt_tot=prhnewprice
replace price_alt_tot=(prhnewprice*1.01) if car_option==`x'
gen price_alt=price_alt_tot/1000

drop price_alt_tot
*/

*Simulate a price increase according to 1% of median price:
gen price_alt=price
replace price_alt=price+0.4644458 if car_option==`x'	

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}

sort id car_option


*Add the mata program here: 

mata {
X=. ;
st_view(X,.,"price_alt height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4");


l_fit=exp(X*beta) ;

/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;

/*
mata drop X1 ncol V
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;
/*
mata drop l_fit X 
*/
l_fit=. ;
X=. ;

prob=mean(prob_all') ;

/*
mata drop prob_all
*/
prob_all=. ;
st_addvar("double", "prob_alt") ;
st_store(.,"prob_alt",prob') ;


prob=.;
}
	
gen d_prob=prob_alt-prob
gen elast=d_prob/prob

frame put d_prob elast car_option, into(elast`x')
frame change elast`x'

gcollapse (mean) d_prob elast, by(car_option)	

gen d_option=`x'
rename car_option s_option

frame change results_elast
fframeappend _all, using(elast`x') drop
*mkmat d_prob, matrix(B`x')

*matrix semi_elast_total[`x',1]=B`x''


frame change default

frame put d_prob elast car_option wealth_group, into(elast`x'_w)
frame change elast`x'_w

gcollapse (mean) d_prob elast, by(car_option wealth_group)	

gen d_option=`x'
rename car_option s_option

frame change results_elast_w
fframeappend _all , using(elast`x'_w) drop


frame change default

	
drop prob_alt price_alt price_w_alt2 price_w_alt3 price_w_alt4 d_prob elast
timer off 1
timer list 1
}	
	
*Now check the matrices and average some of the results: 

*We need the drivetype for each option to make sure we summarise based on the right 
frame change default
frame put car_option drivetype prob choice price, into(car_fuel)

frame change car_fuel

gcollapse (firstnm) drivetype (mean) prob price (sum) car_reg=choice, by(car_option)

frame change default
frame put car_option drivetype prob wealth_group choice price, into(car_fuel_w)

frame change car_fuel_w


gcollapse (firstnm) drivetype  (mean) prob  price (sum) car_reg=choice, by(car_option wealth_group)


	
frame change results_elast


frlink m:1 d_option, frame(car_fuel car_option) gen(d_opt)
frlink m:1 s_option, frame(car_fuel car_option) gen(s_opt)

frget d_fuel=drivetype d_price=price d_car_reg=car_reg init_prob=prob, from(d_opt)
frget s_fuel=drivetype s_car_reg=car_reg, from(s_opt)

gen own=(d_option==s_option)

egen elast_type=group(d_fuel s_fuel), label(elast_type, replace)

drop d_opt s_opt

frame change results_elast_w

frlink m:1 d_option wealth_group, frame(car_fuel_w car_option wealth_group) gen(d_opt)
frlink m:1 s_option wealth_group, frame(car_fuel_w car_option wealth_group) gen(s_opt)

frget d_fuel=drivetype d_price=price d_car_reg=car_reg init_prob=prob, from(d_opt)
frget s_fuel=drivetype s_car_reg=car_reg, from(s_opt)

gen own=(d_option==s_option)

egen elast_type=group(d_fuel s_fuel), label(elast_type, replace)

drop d_opt s_opt


**Save the datasets as results: 

frame change results_elast
save "$root/Results/estimation_results/elasticities/elast_tot_abschange.dta", replace

frame change results_elast_w
save "$root/Results/estimation_results/elasticities/elast_w_abschange.dta", replace


